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Abstract 

We have developed a finite-element micromagnetic simulation code based on the FEniCS 
package called magnum. fe. Here we describe the numerical methods that are applied as well 
as their implementation with FEniCS. We apply a transformation method for the solution 
of the demagnetization-field problem. A semi-implicit weak formulation is used for the in- 
tegration of the Landau-Lifshitz-Gilbert equation. Numerical experiments show the validity 
of simulation results, magnum. fe is open source and well documented. The broad feature 
range of the FEniCS package makes magnum. fe a good choice for the implementation of 
novel micromagnetic finite-clement algorithms. 
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1 Introduction 



Micromagnetic simulations are an important tool for the computational investigation of ferro- 
magnetic materials. In recent years they were successfully used to describe magnetic effects 
ranging from permanent magnets to soft magnetic logic devices to magnetic recording stuctures 
[U 121 |3l H] . Many methods have been proposed to solve the micromagnetic equations numeri- 
cally. Two popular approaches are the finite-difference method combined with the fast Fourier 
transform (FFT) [3 El [7] and the finite-element method combined with the boundary-element 
method [SI El [TU]. Furthermore fast multipole methods |11[ I12j . nonuniform FFT [13] and 
low-rank tensor methods |14[ [T5] have been successfully applied to micromagnetic problems. 

Different open-source codes for the finite-difference method [161 [T71 IT8] as well as for the 
finite-element method |19| [20] are available. Moreover there are a couple of reports on closed 
source simulation tools [21 EU [22]. We present the open-source finite-element code magnum.fe 
that heavily relies on the recently published finite-element software FEniCS |23j and that solves 
the dynamic micromagnetic equations with a combination of two weak formulations. 

Technically challenging tasks appearing in finite-element computations such as numbering 
of degrees of freedom, local to global mapping of cell integrals and numerical integration over 
tetrahedra are handled by FEniCS which offers a variety of finite-element bases including arbi- 
trary order Lagrange elements and produces high performance code for the assembly of system 
matrices. 

The high level of abstraction of FEniCS leads to a very concise code that naturally reflects 
the underlying numerical algorithms. This makes magnum.fe an ideal platform for the im- 
plementation of new micromagnetic finite-element algorithms. Implementing alternative weak 
formulations for certain subproblems can often be done with a few lines of Python. 

This paper is structured as follows. In Sec. [2] we briefly present the theory of dynamical 
micromagnetism. In Sec. [3] we describe the numerical methods that are implemented in mag- 
num.fe, namely a transformation method for the computation of the demagnetization field and 
a weak formulation for the integration of the Landau-Lifshitz-Gilbert equation as proposed in 
|24j . Section [4] gives an overview over the implementation of the algorithms and in Sec. [5] we 
show the validity of our code by means of numerical experiments. 



2 Micromagnetism 

Magnetization dynamics in the framework of micromagnetism are described by the Landau- 
Lifshitz-Gilbert equation (LLGE) 

d t m = — 7(m x H e g) + a(m x dttn) (I) 

where 7 is the gyromagnetic ratio and a > is a phenomenological damping constant that 
depends on the material. The magnetization field m is defined on a domain f2 and assumed to 
be normalized everywhere 

\m(x)\ = 1 , x G O. (2) 

This property is obviously preserved by the LLGE ([I]). The effective field H e g is given by the 
negative variational derivative of the Gibbs free energy U (m) with respect to the magnetization 
m 

(3) 

ett fi M s dm V ' 
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The total effective field H e g is the sum of multiple contributions to the Gibbs free energy 

H c g — H cx + -H"demag -f^zeeman (4) 

where H ex is the exchange field, -Hdemag is the demagnetization field and iif zeeman is a constant 
external Zeeman field. 

The exchange field H ex models the quantummechanical effect of the exchange interaction 
and is given by 

H ex = ^rAm (5) 

where A ex is the exchange constant and M s is the saturation magnetization. Including the 
exchange field in the LLGE ([I]) gives rise to a boundary condition posed on the magnetization 
m 

d v m{x) = , x G dn (6) 

where d u is the normal derivative. This boundary condition is often referred to as Brown 
condition [25]. 

The demagnetization field -Hdemag accounts for the magnetic dipole-dipole interaction. In 
the absence of electric current, the demagnetization field is curl- free and hence can be expressed 
as the gradient of a scalar potential u 

Hdemag = ~ Vll. (7) 

The magnetic scalar potential u itself is the solution of a Poisson problem 

An = M s V • m. (8) 

Boundary conditions for this Poisson problem are given as zero at infinity which is often referred 
to as open boundary conditions. With the Green's function of the Laplace operator the solution 
to ^ can be written in the integral form 

m{x') ■ -, yr^dx' (9) 



u{x = — . j 

47T J \X - X'\ 3 

which directly fullfills the open boundary conditions, see [6]. 

3 Numerics 

3.1 Demagnetization Field 

On a finite region the Poisson equation ^ with Dirichlet boundary conditions is solved by 
the weak formulation 



Vu-Vvdx = M s m-Vvdx V v G V C H 1 . (10) 

^ ^sample 

where the boundary conditions are embedded in the function space V. However the problem 
under consideration has open boundary conditions. Consequently we would have to carry out 



the integration in ( 10 ) over the whole space which is not possible with the finite-element method. 

In order to avoid this restriction we use a method called parallelepipedic shell transformation 
|26j . A finite cuboid shell fi s heii is mapped onto the infinite exterior region f2 ex t G K 3 \f2 cu b kj 



via a bijective transformation X(x) 

X : f^shell — > ^ext- (11) 
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Figure 1: The transformation applied for the demagnetization-field computation, (a) The 
regions in transformed and untransformed space. |(b)| Sketch of the transformation showing the 
construction of the transformation direction. The transformation rays and origins marked by 
gray lines and blue points illustrate the construction of the transformation direction for different 
points in the shell. 



The integration over O ex t = -XX^sheii) arid f2 s heii are connected via substitution by 



f(x(X))dX 



f(x)\det(DX(x))\dx. 



(12) 



When applying this to the left-hand side of (10) the gradients would still be calculated with 



respect to the shell coordinates x instead of the exterior coordinates X. The gradient with 
respect to exterior coordinates X is given by 



[Vxv}(x) = [J~ l V x v](x) 



(13) 



with J = DX{x) being the Jacobian of the transformation. Thus the weak formulation with 
shell transformation reads 



/ Vu-Vv dx + (Vu) T fif Vvdx = M s m-Vvdx V v G V 

'-'^cuboid '-'^shcll '-'^sample 



with the metric tensor g given by 



J _1 ' T J _1 |detJ| . 



(14) 



(15) 



The Dirichlet boundary condition u = on 3(J7 cu b o idU0 s h e ii) is embedded in the function space 
V. The metric tensor g is symmetric positive definite, hence the symmetric bilinear form on 



the left-hand side of ( 14 ) is also positive definite. Thus, by the right choice of the subspace 



V C H 1 , problem (14) has a unique solution. 



3.1.1 Choice of Transformation 

We choose the shape of the shell to be cuboid as described in [26]. The transformation from 



^sheii to Oext is carried out along rays as sketched in Fig. lb In the simple case of a cubic sample 
the transformation has a fixed origin and is radial. However in the general case the origin has 
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to be variable in order to obtain a continuous transformation across shell-patch borders. The 
transformation origin in this case moves on the middle plane that is perpendicular to the shortest 
edge of the cuboid. 

In the following the properties of the one-dimensional transformation in the radial directions 
are discussed. Obviously there are many possible choices for such one-dimensional transforma- 



tions that are bijective and fullfill (11). A suitable transformation distorts the basis functions 



used for discretization in a way that the decay of the potential u may be approximated accu- 
rately. As can be seen in ([9]) the decay of the potential is u oc 1 / 1 as | 2 in the far-field approxima- 
tion. In the following piecewise linear basis functions 



(x) = a + bx 



(16) 



are considered in the untransformed space. In order to obtain a test function decaying with 
1/X 2 the transformation X : O s h e ii — > ^ext has to fullfill 



a! + b 



x 2 



and thus 



X 



x — a 



(17) 



(18) 



Furthermore the transformation X has to map the finite interval [R\, R2] to the infinite interval 



[Ri, 00], see Fig. lb 



X(Rt) 
l/X(R 2 ) 



V 



Ri 



This immediatly yields 



X = R, 



R2 — Ri 



Ri 




R2 - x 

as suitable transformation for linear basis functions. 

When using higher order polynomials as basis functions we use the transformation 

X = R 1 E ^. 

Ri — x 

instead of (21). This transforms second and third order polynomials like 

, 1 



(19) 
(20) 

(21) 



a + bx + cx 



a + bx + cx + dx s ->• a' + b' — + 



/ 1 , / 1 



d 



, 1 



X ' ~ X 2 X s 

which enables a much better approximation of the decaying magnetic potential. 



(22) 

(23) 
(24) 



3.2 Landau-Lifshitz-Gilbert Equation 

We solve the Landau-Lifshitz-Gilbert equation ([!]) with a weak formulation originally proposed 
by Alouges in |24| . In a first step we set v = dtm and multiply with vector test functions which 
yields 



(v — am x v) • <fi dx + / 7(771 x H c s) • 4> dx = V G V 3 
Jn Jn 



(25) 
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The terms of the right-hand side of the LLGE ([!]) and thus also the left-hand side are per- 
pendicular to the magnetization m. Therefore it is sufficient to test the equation with test 
functions (j> £ T m and restrict the solution space of v to T m where T m is the tangent space to 
the magnetization m. Following Alouges we set </> = m x w in (25 ) and restrict the new test 
functions w to the tangent space T m , which yields 



/ (av + m x v) ■ w dx — 7 / H c g(m) ■ w dx = V w G T n 
Jn Jn 



(26) 



This scheme can be extended to an implicit ^-scheme by replacing m with m + 6kv with the 
timestep k and 9 £ [0, 1] 



[av + m x v) ■ w dx — 7 / H c g(m + 



w dx = V w eT„ 



(27) 



In contrast to (25) effective field terms which are linear in the magnetization can be integrated 



implicitly without breaking the linearity of the scheme. E.g. considering only the exchange field 
([5]) and performing integration by parts yields 

2A cx7 



f 2A CX 7 f 

/ (av + m x v) ■ w dx H ^— / V(m + 9kv) ■ Vtu dx = 

Jn VoM s J n 



V w £ T„ 



(28) 



where the boundary condition ^ has been taken into account. When using the whole effective 
field Q it is sufficient to treat H ex implicitly in order to get a stable scheme |13j . This is of 
special practical importance since the calculation of the demagnetization field .ffdemag is very 
time-consuming since the corresponding discretized operators are dense and thus in general not 
feasible to be computed. We therefore calculate the demagnetization potential u as described in 
Sec. 3.1 and treat it only explicitly. Including H cx implicitly and -Hdcmag and H zcem£Ln explicitly 



yields the weak formulation 



ff , \ 1 2A ex7 f 

/ (av + m x v) • w dx -\ / 

Jn W)M S J n 



9kVv ■ Vw dx 



7 / (H z 
Jn 



Vw) • w dx 



2Aex7 

u M 



Vm • Vw dx V w eT„ 



s Jn 



which can be written as 



a(v,w) = L(w) V w E T n 



(29) 



(30) 



where a(v,w) is a bilinear form and L(w) is a linear form. As shown in [24J the bilinear form 
a(v, w) can be written as the sum of a skew-symmetric form and a symmetric positive definite 



form. Thus the problem (29) possesses a unique solution v. 



In order to relieve the tangent-space constraint on the test functions w, (29) has to be 



supplemented by a term that accounts for the part of w parallel to the magnetization m. 
Together with the tangent-space constraint for the solution v the system then reads 



a(v,w)+ / A m ■ w dx = Li 
Jn 



w 



v ■ m a dx = 



V w GV 

V a e V. 



(31) 
(32) 



If (v,X) is a solution to this system then v solves (29). We discretize (31) and the constraint 



(|32j) choosing the same order of finite elements for the scalar field A as for the solution v. This 
leads to a saddle-point problem of the form 



B T O A 



(33) 
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(a) 





Figure 2: Mesh of a spherical sample including the cuboid shell for the demagnetization-field 



computation generated with magnum. fe using Gmsh. (a) Cut through the middle plane of the 



mesh, (b) Explosion view of the the meshed sphere including the cuboid shell patches. 



where A £ ^3^x3^ corresponds to the bilinear form a{v, w), f € M. 3N corresponds to the linear 
form L(w) and B 6 jj 3 -/Vx7V an( j correspond to the integrals in (31) and (32). Since B has 



full rank and A is regular, the system (33) has a unique solution 



A single integration step is carried out by first solving (33) for a given magnetization m(t) 
and then proceed in time by 



rrii(t + k) 



rrii(t) + kvi 
\m,i(t) + kvi 



(34) 



where and Vi are the nodal values at node i of the discretized magnetization m and solution 
v. 



4 Implementation 

Finite element software basically has to solve three sub problems: Mesh generation, system- 
matrix assembly and solution of the resulting linear systems of equations. We use Gmsh [27] 
for mesh generation and FEniCS [23J for matrix assembly and the solution of linear systems. 

For the automated generation of a suitable mesh for the demagnetization-field computation 
we use the C++ interface of Gmsh. magnum. fe is able to produce regular meshes for rectangular 
samples or alternatively read mesh information from a mesh file, which may implement any 
format supported by Gmsh. Then a cuboid shell is wrapped around the sample and meshed 
with Gmsh, see Fig. [2] The number of shell layers which largely influences the quality of the 
demagnetization-field approximation is configurable. 

The system-matrix assembly is done with FEniCS which offers a C++ interface as well 
as a Python interface. For magnum. fe we mainly use the Python interface which leads to a 
very concise code that represents the mathematical problem at hand very naturally. Listing [l] 
and [2] show excerpts of the magnum. fe code, namely the definition of the weak formulations 



Jl4| ) and (|3Tj)-{[32J) with the unified form language (UFL)[28J defined by FEniCS. From these 
form definitions FEniCS creates and compiles a fast C++ code for the matrix assembly via the 
FEniCS form compiler (FFC)[29j. Thus a high performance is achieved although the actual 
programming is done in the scripting language Python. 

However when hitting the limits of FEniCS it is often not possible to extend the functionality 
with Python without facing performance issues. All performance relevant extensions to FEniCS 
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1 # Set up mixed function space 

2 VV = VectorFunctionSpaceCmeshj "CG", 1) # SoLution 

3 VS = FunctionSpace(meshj "CG" , 1) # Lagrange Multipliers 

4 V = W * VS 

5 

6 # Set up test and triaL functions 

7 (v, lamb) = TrialFunctions(V) 

8 (w, sigma) = TestFunctions(V) 

9 

10 # Prefactor for exchange field 

n f_ex = (- 2.0 * Aex * gamma) / (mu0 * ms) 

12 

13 # Set up bilinear form 

14 a = alpha * dot(v, w) * dx \ 

15 + dot ( cross (m, v), w) * dx \ 

16 + sigma * inner(m, v) * dx \ 

17 + lamb * inner(m, w) * dx \ 

18 - 0.5 * dt * f_ex * Dxtvti],]) * Dx(w[i],j) * dx 

19 

20 # Set up linear form 

21 L = f_ex * Dx(m[i]jj) * Dx(w[i],j) * dx 



Listing 1: Weak form of the Landau-Lifshitz-Gilbert equation as explained in Sec. 3.2 A mixed 
function space is used to include the solution v and the scalar field A (called lamb here, since 
lambda is a Python keyword). 



1 # Set up the bilinear form for the demagnetization problem 

2 a = Dx(v, i) * Dx(u, i) * dx(0) + \ 

3 Dx(v, i) * Dx(u, i) * dx(l) + \ 

4 Dx(v, i) * gx[i,j] * Dx(u, j) * dx(2) + \ 

5 Dx( Vj i) * gy[i,j] * Dx( Uj j) * dx(3) + \ 

6 Dx( Vj i) * gz[i,j] * Dx(u, j) * dx(4) 
7 

8 # Create Dirichlet boundary conditions 

9 be = DirichletBC(VSj Constant(9.0) J DomainBoundary( ) ) 

10 

n # Use the SymmetricAssembler to assemble the system matrix 
12 A, An = symmetric_assemble(a , be) 

Listing 2: Weak form of the demagnetization-field problem. The matrix-valued expressions gx, 
gy and gz are multiplied using indices i, j and k and integrated over different subdomains of 
the mesh. 



S 



are thus written in C++. Like FEniCS we use SWIG |30| to exploit the interface of the 
extensions in Python. These extensions include the aforementioned Gmsh interface for the 
generation of the cuboid shells. Furthermore an extension for the assembly of matrices and 
vectors that arise from pointwise calculations of functions was written. This extension works for 
nth-order Lagrange elements and performs calculations on the (auxiliary-) nodes of the function 
space. It is used for the renormalization step (34) which is applied only on the nodes as well as 
for an alternative implementation of the extended system (31)-(32) where the constraint ^ is 
restricted to the nodes. 

The algorithms presented in Sec. [3] are implemented with Lagrange functions, which are 
piecewise polynomial and globally continuous. In case of the demagnetization field the order of 
the elements is configurable. For the solution of the Landau-Lifshitz- Gilbert equation we choose 
1st order elements. Assembly code generated by FEniCS uses Gauss quadrature for integration. 
The degree of quadrature is chosen according to the polynomial degree of the integrand and is 
thus exact for cell-wise polynomial functions. Analytical expressions are interpolated to a given 
degree before integration, e.g. we chose the metric tensor g to be integrated 5th order. 

For the solution of the resulting linear systems of equations FEniCS offers interfaces to a 
variety of open source linear algebra backends. magnum. fe uses the Trilinos Epetra [3T] backend 
for both the demagnetization problem and the solution of the LLGE. Since in both cases the 
resulting system matrices are sparse, iterative Krylov-space methods are applied. The matrix 
of the demagnetization problem is symmetric and positive definite as shown in Sec. |3.1| Thus 
a conjugate gradient solver in combination with an algebraic multigrid preconditioner is used. 
The saddle-point problem arising from the Landau-Lifshitz-Gilbert equation is solved by an 
ILU preconditioned GMRES solver. 



5 Numerical Experiments 
5.1 Demagnetization Field 

For validation and comparison of the demagnetization-field algorithm of different order as pre- 



sented in Sec. 3.1 we choose a homogeneously magnetized unit cube. The energy of this system 



can be computed as 



E = ~ I Vu-mdx (35) 



n 



which is 1/6 if m = (0,0,1). Figure [3] shows results for different polynomial degree of the 
basis functions and different number of non-zero system-matrix entries. The latter was chosen 
as measure since a matrix-vector multiplication is linear in this size and so is the iterative 
solution of the associated linear system. Also the storage requirements are linear in the number 
of non-zero matrix entries. 

Figure [3a] shows the results of the energy calculations. The 2nd and 3rd order method 
perform clearly better than the 1st order method for the same number of matrix entries. This 
is a consequence of the additional 1 /X term that the higher order elements provide in contrast 



to the 1st order elements, see d23| and (24) 



Figure 3b shows the number of iterations needed for the iterative solution of the linear system 
of equations. Together with the numerical complexity of a single matrix- vector multiplication, 
the number of iterations gives the over-all complexity of the demagnetization-field algorithm. 
The log-log plot yields a linear dependence with slope ~ 0.4, resulting in an over-all complexity 
of 0(N 1A ) for the demagnetization-field algorithm. 



Finally Fig. 3c and Tab. 3d show the convergence rates of the algorithm. In order to account 



for auxiliary nodes for the higher order methods the discretization h is not taken from the mesh, 
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Figure 3: Demagnetization energy calculations for different order polynomial test and trial 



functions and different numbers of matrix entries, (a) Energy calculations for a homogeneously 
magnetized unit cube. |(b)| Iterations of the Conjugate Gradient solver, (c) and | (d) | Convergence 
rates for different order polynomials. 
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Figure 4: In-plane magnetization configuration of a so-called s-state in a rectangular thin film 
of the size 500 x 125 x 3nm 3 with the material parameters of permalloy. 

but estimated by iV -1 / 3 where N is the number of degrees of freedom, i.e. the total number of 
nodes and auxiliary nodes. Again the 1st order method shows the poorest performance. The 
2nd and 3rd order methods have convergence rates of approximately 2 which is expected from 
other demagnetization methods, see [7]. 

For further numerical experiments with the presented demagnetization-field method and 
comparsion to other recently developed methods see fTJ. 

5.2 Landau-Lifshitz Equation 

The method for the integration of the Landau-Lifshitz-Gilbert equation is validated with the 
standard problem #4 proposed by the Micromagnetic Modeling Activity Group ^MAG [32] . A 
rectangular sample of the size 500 x 125 x 3nm 3 is relaxed in a so-called s-state with the bulk 
magnetization pointing in the x-direction, see Fig. |4j The material parameters of the sample 
are chosen similar to those of Permalloy 

A cx = 1.3 • 10" n J/m (36) 
M s = 8.0 • 10 5 A/m (37) 
a = 0.2. (38) 

Then, in addition to the exchange field and the demagnetization field, a homogeneous external 
Zeeman field H zeema _ n = (—24.6, 4.3, 0)mT is applied which results in a switching process. Fig- 
ure [5] shows the evolution of the averaged magnetization components in time as calculated by 
magnum. fe compared to the results of the finite-difference code MicroMagnum |18j . 

6 Conclusion and Outlook 

We present the open-source micromagnetic software magnum. fe. magnum. fe is a complete three- 
dimensional finite-element code, which computes magnetization dynamics with a combination 
of two linear weak formulations. It is written in C++ and Python and uses the finite-element 
package FEniCS [23J . The correctness of the code is demonstrated by a number of numerical 
experiments including the fj,MAG standard problem #4. 

Due to the multitude of features and the high level of abstraction of FEniCS, magnum. fe is 
well suited for the implementation of novel finite-element algorithms, magnum. fe itself is well 
documented and unit tested and is freely available at github [33J. 

We plan to implement alternative demagnetization-field algorithms as well as integration 
schemes for the Landau-Lifshitz-Gilbert equation. Contributions to magnum. fe are very wel- 
come. 
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Figure 5: Time evolution of the averaged magnetization for the standard problem $4. The 
results of magnum. fe are compared to that of the finite-difference package MicroMagnum. The 
different branches correspond to the components of the averaged magnetization. 
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